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I. INTRODUCTION 



Let T be the life time of a component that is 
exponentially distributed with parameter A. Because of 
the manufacturing process and different logistical forces 
acting on the components, A behaves as a random variable 
with a certain distribution G. The distribution of the 
time to failure, under the above assumptions, is called a 
mixed exponential distribution. 

Let N components be put on test until a test time T Q 
is reached. Suppose that as a result of testing, 
tijtg > ♦ • • jt^ were life times of components that had failed 
and t k+1 ,t^ + 2 , . • . ,t N were life times of components that 
had not failed. If G is a Gamma distribution with shape 
parameter a and scale parameter B, then a and B can be 
estimated by using the method of maximum likelihood. 

The M.L.E. of the reliability function of the components 
can be expressed as 



R(t) = ^ — x- 

(1 + Bt ) a 



A 

where a is the M.L.E. of a 

A 

B is the M.L.E. of B 



Computer simulations were conducted to determine the 



sensitivity of this estimation procedure when the exponential 



distributions are mixed differently from the assumed Gamma 



distribution mixture. 
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II. STATEMENT OP THE MODEL 



The model and the derivations In this section are 
discussed in detail in a paper by Myhre and Saunders [ 5 ]. 

Let T the life time of a component be exponentially 
distributed with parameter A, where A is a random variable 
with Gamma distribution G with shape parameter a ana scale 
parameter 3. 

The density function of A is 



g(X) 



dG(A) 



A a -V X/g 

r(a)e a 



(1) 



The reliability of the component is 



R(t) = P (T>t ) 



= / P (T>t | A= A ) dG ( A) 
0 



/ 

0 



e 



— At 



r (a) 3 a 



dA 



1 

(1 + 3t) a 



f 

0 






a-1 

A 



r(a) ( 



1+Bt 



) a 



dA 



1 

(1 + 3t) a 



( 2 ) 
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Then the density function of T is given by 



f(t) = - 



dR(t ) 
dt 



( 1+Bt ) a+1 



( 3 ) 



Suppose that of N items that were put on test it was 
observed that t^ , tg , . . . , t^. were the life times of k items 
that had failed before T Q and let t k+1 = = ^o 

were the life times of items that had not failed. The 
maximum likelihood estimators of the parameters are derived 
as follows: 

The likelihood function is 



k 

L(ct , 6 ) = n f(t.) 

i=l 1 



j =k+l 



N 

n r( t . ) 

V4-1 J 



(*n 



Substituting equations (2) and (3) into (A) 




j=k+l (1+Bt . ) a 

J 




k 

n d+Bt.) 

i=l 1 



N 

n (i+Bt . ) a 



j =k+l 




(5) 



N 



n (i+Bt . ) n (i+Bt,) a 
i=l 1 i=l 1 
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Taking the logarithm of L(a,8), 



L n = In L(a,B) 



N 

= k ln(a8) - a Z ln(l+Bt.) - Z 

i=l 1 i 



The maximum likelihood estimators of a and 
obtained by differentiating L n 



8a 



k 

a 



N 

Z ln(l+8t.) 
i=l 1 




k 

8 



a 



N 

Z 

i=l 



l+8t ± 



k 

Z 

i=l 



1+Bt 



i 



A 

Suppose first that a is known. Then 8, 
of 8 can be computed by setting equation (8) 



k 

8 



N 

a Z 
i=l 



t . 

i 



l+8t ± 



k t A 
i=i 



0 



ln( 1+Bt . ) 

1 1 

( 6 ) 

can be 



( 7 ) 

( 8 ) 

the estimate 
to zero 

( 9 ) 
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Rearranging equation (9), 



k N t ± k t ± 

6 a i f 1 i+et 1 + 1 f 1 i+Bt i 



N Bt 

k = a Z ?Tor - 
i=l 1+6t i 



Bt. 



+ Z 



i=l 1+6t i 



N 

= a Z (1- 
i=l 



1+Bt 1 ) + i l 1 (1 1+Bt 1 ) 



N 



= Na - Z 



i-1 1+Bt i 



k 

+ k - £ 



i=l 1+Bt i 



Finally , 



N 

a Z 



k 

+ Z 



i=l 1+Bt i i-i 1+Bt i 



- Na = 0 



( 10 ) 



The solution of equation (10) is B, an estimate of B. 
Note that since 



N 



A(B) = a Z 



+ Z 



1=1 1+6t i ' k^i 1+6t i 



- Na 



is a decreasing function of B with A(0) = k and 

lim A(B) = -Na, there exists exactly one solution for B 

0-*co 

whenever a > 0 and N >_ k >_ 1. 

A 

If B is known and B is positive then a could be computed 
directly by setting equation (7) to zero, 



Z ln(l+Bt . ) 
i=l 1 



(ID 
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If a and 8 were both unknown, a and 8 could be obtained 
by solving the equations simultaneously 



k 

a 



N 

l ln(l+8t . ) = 0 
i=l 



( 12 ) 



t r - a 



N t ± 
1=1 1+Bt i 



k 

- E 



i-1 1+6t i 



= 0 



(13) 



Eliminating a from equations (12) and (13) results in an 
equation that contains only 8, say B(8) = 0 where 



1 N k 

B ( 8 ) = £ E ln(l+8t . ) Z v+bV 
k i=1 i i=1 1+BtjL 



N 

- E 



8t. 



i=l 1+et i 



(14) 



Then 8 is the solution of B(g) = 0 and a can be solved 
using equation (11). 

We now proceed to examine the behavior of the function 

A 

B(8) . In particular 8 exists only under certain conditions. 
Note that B(0) = 0 and lim B(8) = -N and B ( 8 ) is a continuous 

g-voo 

function. Hence, B(8) = 0 would have a positive solution 
if its derivative at 8 = 0 is positive. Since B’(0) = 0, 
the sign of lim — £ - - may be examined. A condition that 



8-*-0 td i / p \ 

B(8) = 0 has a positive solution is that lim — p v ■ > 0 

8-^0 B 

(see Figure 1). 



8 



11 



Dl 

m 

mc 



-Al 



Figure 1 
Graph of B(B) 



Rewriting equation (14) 



, N k .. N Bt. 

b(b) - £ z m(i+Bt ) z Ti ± F - - z TT ^- 

k 1=1 i 1=1 1+Bt i ±=1 l + pt ± 

1 ^ ^11 ^ 

B ’ (e) = e ^ 1^7 i4t7 - r ln(1+6t i> ^ 



N 

- Z 



N 

+ Z 



Bt j _ 



N 

= Z 



i=i 1+Bt i i=l (l+Bt 1 ) 2 
2 



Bt, fc , N 

- r- Z 



i=l (1+Bt ± ) 2 k i=l 1+8t i 



t k l 

1 [k - Z X 



i=l 1+Bt 



, N k t. 

- £ Z ln(l+3t.) Z i — p 

K i=l i=l (1+Bt i ) 



(1+Bt.) 2 



- ] 
i 

(15) 
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Bill) = 

B 



N t . 

E , 

i=l (l+3t i )' 



1 N t^ k 

k 1 f 1 l+6t i i f 1 l+6t 1 



, N ( 1+Bt. ) k t, 

- £ E in g-i- E i— 

1=1 p 1=1 ( l+6t ± ) 



( 16 ) 



Note that lim 
6+0 

( 1 6 ) we get 



ln(l+Bt ± ) 

B 



t^. Taking the limit of equation 



pi f o\ N p pN k 

lim -- a - 1 = E t/ - r E t . E t, . (17) 

B +0 6 1=1 1 k 1=1 1 1=1 1 

B ? ( g ) 

Since it is required that lim — 1 — - > 0 to get a positive 
- B+0 15 

solution for B, it can be stated from equation (17) that 

A 

a sufficient condition for the existence of a positive B 
is that 



N p N k 

E t. > £ E t. It. (18) 

1=1 1 K i=l 1 1=1 1 



The reliability estimate R(t) could be computed using 

A 

equation (2) whenever a positive solution for B exists. 

A A 

If a positive solution for B does not exist however, R(t) 
can not be computed using this equation. In this latter 
case, the reliability could be estimated in the following 



manner. 



Assuming A is Gamma G(a,3) distributed, then E(A) = aB. 
Setting aB = Y and substituting a = y/B into equation (6), 
the log likelihood function becomes 
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( 19 ) 



N k 

L s k In y - j Z ln(l+Bt.) - Z ln(l+3t,) 
n P i=l 1 1=1 1 



Differentiating L n we get 



3L n 

3Y 




ln(l+Bt 1 ) 



( 20 ) 



9L n 



X 

B 



N 

Z 

1=1 





N 

I 

1=1 



ln(l+Bt i ) 



l=l 1+Pt l 



. ( 21 ) 



The M.L.E. of y and B could be obtained by setting equations 
(20) and (21) to zero, and solving simultaneously for y and 
3. The resulting estimator for y is 




Z ln(l+3t,) 
1=1 



Substituting y from equation (22) into equation (21) and 
rearranging , 



1 

k 



k 

Z 

i=l 



1 

1+Bt ± 



N 

Z 

i=l 



ln(l+Bt 1 ) = 



N 

Z 

i=l 



Bt ± 

1+BtT 



which is exactly the same as equation (14). Thus the M.L.E. 
of 3 is the same as before. However whenever a positive 

A 

solution for 3 does not exist the reliability estimate 
could be computed by treating the life times of the com- 
ponents as exponential random variables with a constant 



01 

Hi 

Mi 



failure rate y, that could be computed by taking the 
limit of equation (22) to get 



Y 



lim 

6-0 



k0 

u 

E ln(l+0t . ) 
i=l 1 




In this case, the reliability estimate is given by 
R(t) = e"^ 



( 23 ) 



( 24 ) 
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III. COMPUTER SIMULATION 



To determine the robustness of this estimation procedure, 

computer simulations were conducted for a discrete probability 

distribution on A. Let A_,A-,...,A be the chosen values 

12 m 

of failure rates of the components with associated proba- 
bilities p^ ,? 2 » • • • >P m * A total of N exponential failure 
times T's were generated using a standard exponential 
random number generator. The T’s were then truncated at 
a planned test time Tq. If this set of T’s gave any posi- 

A A 

tive solution for B, then B could be computed by setting 

A 

equation (14) to zero and solving for B and a could be 
obtained from equation (11). The reliability estimate was 
computed by using 



H(t) = ( 25 ) 

(1+Bt ) a 



If no positive B exists, then the exponential procedure 

A 

was used, y the estimate of E(A) was computed by using 

Y = -k (26) 

Z t + (N-k)T 
1=1 
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where 



N: the number of components on test 

T q : the planned test time 

k: the number of components that failed before Tq. 

The reliability estimate then was computed using equation 

(2*1). 

Since the true reliability of the component is 
m -A, t 

R(t) = £ p.e 1 (27) 

i=l 1 

where m is the number of chosen failure rates, a measure 

A 

of robustness of this estimate is R(t) - R(t). 

A delineation of the simulation steps is as follows: 

(1) Generate N times to failures T^jTg,...,^ where 
T^ is exponentially distributed with parameter 

X i and each A^ is randomly generated in accordance 
with the stated probability distribution on A. 

(2) Truncate T^,T 2 ,...,T^ at a test time Tq. 

(3) Check if this set of T^,T 2 ,...,T^ give any positive 

A 

3 by taking a small value Bq then check the value 
of B(Bq) from equation (1*0. 

If B(Bq) is positive, continue to step (*0. 

A 

If B(Bq) = 0, then B = Bq, and continue to step (5). 
If B(Bq) < 0, use equation (26) and (2*1) to compute 
the reliability estimate then continue to step (7). 
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(4) Compute 3 from B(3) = 0 using Newton-Raphson method. 

A 

(5) Compute a using equation (11). 

A 

(6) Compute the reliability estimate R(t) using 
equation (25). 

(7) Replicate steps (1) through (6) r times. 

(8) Compute the mean value and standard deviation of 

A ^ 

all values of R(t), and let R(t) be the mean value 

A 

Of R(t ) . 

(9) Compute the true reliability R(t) by using equation 
(27). 

A 

(10) The measure of robustness is R(t) - R(t). 

A second simulation was done as a comparison to the 
first one. The procedure was similar to the first one 
except that instead of truncating the failure times 
T*i ,T^ , . . . ,T^ at a test time Tq (see step 2), they were 
truncated at the time of the failure. In other words 

N components were put on test, and the test was terminated 
after k components had failed. 

A 

Using this procedure however, a positive solution for 3 
will never exist for k £ 2. Only if k > 3 under certain 

A 

conditions will a positive 3 exist. This can be shown as 
follows . 

Recall that a sufficient condition for the existence 

A 

of a positive 3 (see inequality (18)) is. 
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N ? N k 

Z t, d > £ Z t, Z t, 

1=1 1 K 1=1 1 1=1 1 

For k=l then t ^ =t N-l = * * * =t 2 =t 1 an< * the ine Q ualit y ( 1 8 ) 

becomes 

Nt x 2 > 2Nt 1 2 

and this, of course is invalid. 

For k=2 , t^=t^_-^= . . . =t 2 an< ^ the i ne Q ua lity ( 1 8 ) 
becomes 



t ± 2 + (N-l)t 2 2 > t x 2 + (N-l)t 2 2 + Nt x t 2 
and this is also invalid. 

For k=3, t N =t N-l = ’ * ‘ =t 3’ the t ne 9 ual tty ( 1 8 ) becomes 

t l + t 2 2 + ( N " 2 ) t 3 2 > |(t 1 +t 2 +(N-2)t 3 )(t 1 +t 2 +t 3 ) 
or 

^[t 1 2 + t 2 2 + (N-2)t 3 2 ] > |[2t 1 t 2 +(N-l)t 1 t 3 +(N-l)t 2 t 3 3 
and this inequality could be valid for some sets of t^'s. 

A 

Thus, using this' later simulation procedure 6 > 0 will 
possibly exist only for k ^ 3. 

Using these two procedures, simulations were done and 
the results appear in Section IV. 
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IV. RESULTS OP THE SIMULATION 



Graphs which depict the results of the simulation 
appear on the following pages. Symbols used on the graphs 



are : 




R(t): 


The reliability at time t 


R(t): 


The estimate of R(t) 


N: 


Number of components on test 


T : 

o 


The planned test time 


r : 


Number of replications 


A: 


Parameters of generated life time T's of the 
components 


p: 


Discrete probability distribution for A 


k: 


Number of components that failed before T Q the 
planned test time 


PB: 


Number of replications that gave positive 
solutions for 6 


NB: 


Number of replications ^that did not give a 
positive solution for 8 


Total : 


: Total number of replications that had k failures 



To illustrate the graphs, let us observe Figure 2. This 
graph shows that five exponential life times (N=5) were 
generated via a set of A: .010, .010 with probability 

density p: .50, .50. By using a planned test time T q = 10, 

the reliability estimate was computed. With 50 replications 

A 

(r=50) the average reliability R(t) was computed and was 
shown on the graph. Further, it can be observed that for 
k = 0 (no failure occurred), NB = 28 (28 replications 
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occurred without any failure). It happens that 3 replica- 
tions occurred with 2 failures (Total = 3, k = 2) and all 

A 

of these three replications gave no positive 8 (NB=3) and 

A 

of course no replication gave positive solutions for 6 
(PB=0 ) . 

Another example is Figure 25 where 50 replications 
(r=50) of 10 exponential life times (N=10) were generated 
via a set of X: .005, .020 with probability density 

p: .50, .50. It is observed that by using a planned test 

time T q = 20 unit times, 22 replications (Total=22) occurred 
with 2 failures (k=2) of which 5 replications gave positive 

A 

solutions for 8 (PB=5) and the rest did not give any 

A 

positive solution for 8 (NB=17). 

Figure 4*4 shows that 50 replications (r=50) of 
10 exponential life times (N=10) were generated via a 
set of X: .005, .010, .015, .020 with probability density 

p: .*40, .20, .35, .05. The N items were tested until a 

fixed number of failures (k=5) occurred. Eight of those 

A 

replications gave positive solutions for 8 (PB=8) and the 

A 

rest (NB=*42) gave no positive solution for 8. 
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NB: 




2 




7 


15 


6 


3 


2 








Total : 


2 




12 


18 


13 
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2 
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Figure 39 

N = 20 
T 0 = 20 
r = 30 



X: 


.005 


.007 


.010 


.015 


.017 


.020 .025 


.027 


.030 


.035 


p: 


.25 


.20 




.20 


.10 


.07 


.05 


.05 


.03 


.03 


.02 


k: 


0 


1 


2 


3 


4 


5 


6 


7 


8 








PB: 


0 


2 


1 


2 


4 


1 


2 


1 


1 








NB: 


1 


3 


5 


6 


6 10 


2 


2 


1 








Total : 


1 


5 


6 


8 


10 11 


4 


3 


2 
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Figure 40 



N = 50 
T q = 20 
r = 50 



X: 


.005 


.007 


.010 


.015 


.017 


.020 


.025 




C\J 

O 

• 


o 

on 

o 


.035 


p: 


.25 


.20 


.20 


.10 


.07 


.05 


.05 




.03 


.03 


.02 


k: 


4 


5 


6 


7 


8 


9 


io : 


11 12 


13 


1H 


15 






FB: 


1 


0 


1 


1 


0 


3 


l 


1 


0 


0 


1 


1 






NB: 


1 


0 


1 


i| 


4 


2 


7 


3 


6 


5 


6 


1 






Total: 


1 


0 


1 


5 


il 


5 


8 


4 


6 


5 


7 


2 
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Figure 4l 



N = 10 
r * 50 

X: .005 - 02 

p: .50 -50 

k: 3 

PB: 4 

NB: 46 

Total: 50 
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'97 
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Figure 42 

N = 10 
r = 50 

X : .005 -020 

p: .50 .50 

k : 5 

PB: ' 8 

NB: 42 

Total: 50 
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Figure 4 3 



N = 10 
r = 50 



X: .005 .010 .015 .020 

p: .40 .20 .35 .05 



k: 3 
PB: 4 
NB: 46 
Total: 50 
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N = 10 
r = 50 

X : .005 .010 .015 .020 

p: .40 .20 .35 .05 

k:' 5 

PB: 13 

NB: 37 

Total: 50 
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Figure 45 

N = 10 
r = 50 

A: .005 .007 .010 .015 .017 .020 .025 .027 .030 .035 

p: .25 .20 .20 .10 .07 .05 .05 .03 .03 .02 

k: 3 

PB: 4 

NB: 46 

Total :50 
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Figure 46 



N = 10 
r = 50 

X: .005 .007 -010 .015 .017 .020 .025 .027 .030 

p: .25 .20 .20 .10 .07 .05 .05 .03 .03 

k = 5 
PB = 9 
NB = 41 
Total =50 



.035 

.02 
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V. CONCLUSION 



The results of the simulations show that 

1. In all cases, under the first procedure of 

A 

simulation a positive solution for B was never obtained 
for relatively small planned test times, e.g., T Q = 10. 

A 

Some positive solutions for 3 were obtained when T q = 20. 

2. By using the second procedure (fixed number of 

A 

failures) a positive solution for B never existed for 
k = 1 and k = 2. Only for k >_ 3 were positive solutions 

A 

for B obtained. The larger the value of k, the more 
accurate the computed reliability estimate will be. 

3. The result also showed that most of the reliability 
estimates have quite small error relative to the true value 
of the reliability. It was less than 1$ on the average. 

The smallest relative error of the reliability estimates 
was .02$ while the largest one was about 3*15$. Thus the 
proposed estimation procedure is quite accurate when the 
underlying distribution on A is a discrete distribution 
similar to those assumed rather than a continuous gamma 
distribution. In this regard the estimation procedure 
indicates some property of robustness. 
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SYMBOLS USED IN THE COMPUTER OUTPUTS 



N: Number of generated component's life times 

T q : Planned test time 

Parameter: Parameters (X) of the life times 

Probability: Discrete probability density of X 

k : Number of failure components on test 

PB: dumber of replications that gave a solution for 

3 > 0, associated with k 

NB: Number of replications that did not give a solution 

for 3 > 0, associated with k 

Total: Number of replications that gave k failures 

T: Time to compute the reliability and its estimate 

R(T): True value of the reliability 

RE(T): Reliability estimate of R(T) 

SD: Standard deviation of 50 replications of RE(T) 

MOE: Measure of Effectiveness of estimation, i.e., 

RE (T) - R(T) 
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COMPUTER OUTPUT 
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